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Abstract 

The  two-  and  three-  dimensional  development  of  leading-edge  receptivity  and  transition  in 
a  2-D  Joukowsky  airfoil  boundary  layer  are  investigated  by  direct  numerical  simulation  (DNS) 
using  the  incompressible  Navier-Stokes  equations.  The  numerical  investigation  are  based  on  the 
so-called  spatial  approach.  A  contravariant  based  governing  system  is  derived  so  that  we  can 
simulate  the  receptivity  and  transition  with  complex  geometries.  The  numerical  results  agreed 
very  well  with  the  experimental  results  for  the  flat  plate  case,  and  also  agreed  well  with  the 
result  obtained  by  other  researchers  for  the  2-D  elliptic  leading-edge  receptivity  case.  Some  new 
phenomena  for  the  transition  around  Joukowsky  airfoils  were  observed  which  at  least  qualitatively 

agree  with  physics. 
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1  Introduction 

To  understand  the  origins  of  transition  and  turbulence  remains  a  critical  problem  of  fluid 
mechanics.  Important  sources  of  information  regarding  applications  of  transition,  for  example,  the 
design  of  laminar  flow  control  systems  to  delay  or  even  eliminate  boundary  layer  transition  and  to 
increase  the  aerodynamic  performance  of  the  systems,  can  be  found  in  a  variety  of  recent  literatures 
(e.g.,  AGARD  Report  No.786,  1992;  No.793,  1993).  Receptivity,  the  process  through  which 
external  disturbances  generate  instability  waves  in  boundary  layer  (Morkovin,  1969),  provides 
the  vital  initial  conditions  of  amplitude,  phase,  frequency,  and  wavelength  for  the  breakdown  of 
laminar  flow.  Clearly,  the  study  of  receptivity  promises  significant  advances  in  practical  transition 
prediction  methods.  This  has  been  the  subject  at  the  first  stage  of  our  project. 

The  process  of  boundary  layer  transition  can  be  qualitatively  described  as  the  following  sce¬ 
nario:  receptivity  to  freestream  disturbances,  initial  growth  of  instability  waves,  three-dimensional 
and  nonlinear  effects  coming  into  play  and  thus  leading  to  breakdown  to  turbulence  (Saric  &  Reed, 
1994).  As  briefly  reviewed  by  Liu,  Liu  &  McCormick  (1993),  considerable  progress  has  been 
achieved  with  theoretical,  computational,  and  experimental  tools  on  broadening  our  knowledge 
of  the  latter  two  stages  of  transition  for  boundary  layer  in  the  last  two  decades.  However,  efforts 
on  understanding  the  mechanism  of  receptivity  are  still  quite  limited.  Mathematically  speak¬ 
ing,  receptivity  problem  differs  from  the  stability  analysis,  which  describes  the  normal  modes 
(eigenvalues)  of  disturbances  determined  from  the  linearized  N-S  equations,  due  to  the  fact  that 
its  governing  system  being  typically  the  full  N-S  system  with  appropriate  boundary  and  initial 
conditions.  Because  of  the  system’s  inherent  susceptivity  to  many  subtle  environment  influences, 
receptivity  theory  remains  rather  incomplete  with  regarding  to  quantitative  predications  of  the 
strength  of  the  disturbances  introduced  into  the  boundary.  Well-documented  experiment  data 
in  this  area  are  also  limited.  DNSs,  on  the  other  hand,  are  playing  increasingly  important  role 
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in  the  investigation  of  transition  with  the  recent  progress  made  in  the  development  of  extremely 
powerful  computers  and  numerical  algorithms  (Reed,  1993,  and  Kleiser  k  Zang,  1991).  It  is  men- 
tionable  that  significant  achievements  in  spatial  DNS  have  been  made  by  Liu,  Liu  k  McCormick 
(1991,1992,1993),  and  Liu  k  Liu  (1993,1994,1995)  with  the  support  of  NASA  Langley  Research 
Center  for  1990-1994  and  US  AFOSR  for  1994-present.  Our  former  and  present  DNS  researches 
have  established  a  viable  framework  for  the  study  of  boundary  layer  receptivity. 

Among  various  paths  through  which  external  energy  is  internalized  in  form  of  instability  waves 
in  the  boundary  layers,  such  as  curvature,  sweep,  roughness,  freestream  disturbance,  and  surface 
pressure  gradient,  we  choose  to  consider  the  receptivity  to  freestream  vortical  disturbances  of 
boundary  layer  flows  which  include  the  boundary  layer  flow  over  a  semi-infinite  flat  plate  with 
an  elliptic  leading-edge,  and  that  over  a  Joukowsky  airfoil.  This  is  referred  to  as  leading-edge 
receptivity.  Related  works  have  been  performed  by  Lin,  Reed  k  Sane  (1992),  and  Buter  k 
Reed  (1994)  which  showed  the  clear  appearance  of  the  T-S  modes  in  the  boundary  layer  over 
a  two-dimensional  semi-infinite  flat  plate  with  elliptic  leading-edge  of  different  aspect  ratio  and 
some  other  modified  geometries.  Receptivity  to  freestream  spanwise  vorticity  is  found  to  be 
different  for  two  cases  considered:  symmetric  and  asymmetric  forcing  of  the  same  strength  of 
input  disturbances.  For  the  present  work,  the  main  purpose  is  to  validate  the  capability  of 
our  new  spatial  DNS  code  which  is  motivated  for  further  studies  of  real  time-dependent  three- 
dimensional  compressible  turbulent  flows  around  Air  Force  flight  vehicles  (the  preliminary  study 
is  limited  to  incompressible  transitional  flow).  We  first  study  the  receptivity  of  Joukowsky  airfoil 
to  freestream  vortical  disturbances.  So  far,  all  numerical  studies  on  receptivity  are  limited  to 
2-D  simple  geometry,  and  serve  for  the  purpose  of  the  receptivity  study.  This  work  is  a  first  try 
to  study  receptivity  for  3-D  and  more  general  configurations.  Finally,  we  investigated  the  whole 
process  of  flow  transition  around  Joukowsky  airfoils. 


2  Governing  System 


The  3-D  incompressible  time-dependent  Navier-Stokes  equations  are  considered  as  the  govern¬ 
ing  system  in  this  work.  On  the  general  curvilinear  coordinate  system  (f,r?,C),  these  equations 
can  be  expressed  as  follows: 
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Unfortunately,  equations  (l)-(3)  can  cause  some  trouble  for  the  numerical  simulation  when 
the  perpendicularity  of  U-,  V-  and  W-  are  fairly  poor.  In  order- to  give  a  system  with  better 
behavior,  we  reorganize  them  to  get  equations  along  U-,  V-,  and  W-  directions,  and  named 
them  as  U-momentum,  U-momentum,  and  ^-momentum  equations.  There  are  two  ways  to 
obtain  the  discretized  system: 


1.  Obtain  the  U— ,  V  — ,  and  W—  equations  analytically,  then  discretize  them,  or 
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2.  Discretize  equations  (l)-(3)  first,  and  then  combine  them  to  obtain  the  discretized  U-,  V  , 
and  W—  equations. 

In  the  first  way,  a  non-derivative  term  will  be  generated  in  each  momentum  equation  which  makes 
the  conservation  of  system  very  poor,  and  sometimes  it  is  difficult  to  maintain  the  zero  residual 
requirement  of  the  momentum  equations  for  the  uniform  flow  field.  Therefore,  we  choose  the 

second  way. 

To  minimize  the  effects  of  machine’s  round-off  error,  we  decompose  the  total  flow  into  two 
parts:  base  flow  and  perturbation. 


2.1  Basic  state 

Using  the  second-order  central  difference  as  an  example,  we  show  how  to  derive  the  discretized 
governing  system.  First,  we  discretize  (l)-(3)  as  follows  (see  Figure  1) 

WfUf  —  WbUb 


1  dll .  Ueue  UWUW  Vn 
( J~dt',:,k+  AC  At? 


+ 


AC 


1  dv„  .  Ueve  -  Uwvw  ,  Vnvn-Vsvs  ,  WfVf  -  WbVb 


■  + 


AC 


AC  ‘  At? 

+[7(c4+7? 4+c4)F]ijfc_(^AlU),ifc  =  °’ 

1  dw  UeW^-UyU Te  VnWn  ~  .  WjWj-WbWb 

(jUt *tjk+ 

+t7^4  +  7?4  +  c4)F]^-(7iAlU,)^  =  °- 


(9) 


(10) 


(ii) 


Figure  1.  Stencil  in  the  (C,  if)  plane. 
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,  X  MG'! +  X  (11),  and  assuming  u,  v, 

To  derive  the  U-e quation,  using  Cx0*  x  (9)  +  £,ytjk  x  (iuj  ■+■  C2,jfc  V 

defined  everywhere,  we  obtain: 
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(15) 


SviSijk  (<f>) 


-Uw[uw{<f>Xijk  -  4>xw)  +  Vw{4>yijk  ~  <t> yj  +  ww{<Pzi]k  ~  &«)]} 

+  -^—{Vn[un(4>Xi]k  ~  <j>x„)  +  Vn{4> yijk  -  4>yn)  +  wn{4>zi}k  -  4>Zn)] 

Arj  J 

-Vs[us(<t>Xijk  -  <fi*s)  +  v‘(4>y,}k  -  K)  +  MKk  -  ^)]> 

+  -^{Wj[uf{4>Xijk  -  <pXf)  4  Vf(<t>y.jk  -  <t>y,)  +  vMliik  ~  ^ 
-Wb[ub{4>Xijk  -  4>xb)  +  Vb(<t>y,jk  -  <Pyb)  +  wb{4>zijk  -  fab)}}’ 

=  ~^{^l [Ui+lMfeijk  ~  ^*4+lifc)  +  ~  + 

Vi+ljk{4>yt}k  ~  ^Vi+ijk)  +  vi-ljf=(4>yi}k  ~  ^Vi-ijk )  + 

Wi+ljk(<f>zijk  ~  4>zi+xik)  +  wi-ljk(4>z,,k  ~  <t>Zi- ljffc)] 

4  ir\\Uii+^k{4>xi]k  -  <t>x,}+lk)  +  uij-lk(<t>xijk  -  <t>Xij- 1*)  + 
Vij+lk{4>ytjk  -  ^y.j+u )  +  vij-lk{<t>y,]k  -  $ Vij-lk )  + 

Wij+lk{<f>zijk  ~  <t>Zij+ik)  +  wii-lk{<l>zijk  ~  4>Zij-ik )] 

+  |  [Uijfc+l  (<£x,7k  “  <t>x,]k+1 )  +  “tjfc-l  “  ^xUk- 1 )  + 

U«iJfc+l  (^yo*  “  ^V.jfc+l)  +  vijk-l{4>yijk  -  Qytjk-l)  + 

Wijk+l{<f>ztjk  -  4>Zijk+ 1)  +  Wijk-l{<Pz,jk  -  <t>zijk- 1  )]}• 


(16) 


2.2  Perturbation  form 

The  fourth-order  central  difference  formulation  for  the  perturbation  flow  can  be  derived  as 


follows.  Using 

10«  .  3[(tf+t/fc)u +  «*£/]  ,  a[(V  +  Vt)«  +  ^V]  ,  d[{W  +  Wb)u  +  ubW] 
lli+  d£  dr,  d( 

+7(Cl^  +  ^  +  C*I;)F  "  Jr^AiU 

i  ,  a[(u  +  Ui>)^  +  »fc^]  g[(u  +  H)«  +  »bU]  a[(w  +  wb)v  +  vbW] 

7~dt  'r  dr, 

+7«4+4+c4)p-^Ai" 

1  flu,  arcu  +  I7fc)ti»  4  u^t/l  .  5[(P  + Vi)iD  +  tBfcV]  ,  a[(W4  Wfc)n>4  wbW] 

71i  + - T( - +  5v  +  »c 

4«4+,j4+c4)p-^i» 


0,  (17) 


0,  (18) 


0,  (19) 


where  (t*,  vb,  wb)  and  (Ub,  Vb,  Wb)  are  base  flow  velocity  components  and  their  associated 
contravariant  velocity  components,  we  have 


ildu\  ,  1 

[J  dtUjk  12A£ 


[—  (Ui+2jk  4"  Ubt+2jk  )Ui-r2jk  4  ${Ui+ljk  4  Ubi+ljk)ut+ljk 
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-HUi-ijk  +  Ub,_ljk)ui-ijk  +  (U,-2jk  +  Ubi_2)k)Ui-2jk]+ 

— TT-[-(^i+2fc  +  VbiH2k)UH+ +  8(V*'i+lfc  +  VbiJ+lk)Uij+lk 

12At? 

-8(Vij- ik  +  Vbij_lk)uij-ik  +  {Vij-2k  +  Vblj_2k)ui^2k}  + 

— -[-(Wijk+2  +  Wbijk^)uijk+2  +  8(Wii*+l  +  ^+>.^+1 
1 2  ^ 

-8  {Wijk-l  +  WVi)^-!  +  (W'3k~2  +  W^ik-2)ui}k- 2]  + 

-Ui+2jkUbi+ijk  +  SUi+ukUbMih  -  Wi-UkUbi-w  +  U^jkUbi-2^  + 

—  “  '  12A| 

-V-j+2fcUb,.  .+afc  +  Wij+ikUbii+lk  -  + 

12A?7 

-Wtifc+2Mti,-fc4.a  +  %Wtjk+lUb„k+1  ~  8^-1^-.  +  Wijk-2Ubijk_2  + 

"  —  12AC 

{j&T(  +  ^  +  Cl|;)p}"*  -  (7  TeA'u)<it  = 
+  12^[~(^+2jfc  +  UbxJr7}k)Vi+2jk  +  HUi+ljk  +  Ubi+ljk)vt+l]k 
—8(Ui-ijk  +  Ubi_ljk)vi-ijk  +  ( Ui-2jk  +  Ubl_2jk)v,-2jk]  + 

— l—  [-{Vij+2k  4-  Vbij+2k)vij+2k  +  8{Vij+ik  +  VbtJ+Ji\3+ih 

12A7J 

-8  {Vi3-lk  +  Vbi}-lk)Vij-lk  +  (Vi]-2k  +  VbtJ-2k)vi]-2k]  + 

— T~, :[—  {Wijk+2  +  Wbijk+2  )vijk+2  +  8(M'rijfc+i  +  Wbijk+ i)vxjk+i 

— 8(H/ijfc_i  +  W^K-k-1  +  (Wijk-2  +  Wbljk_2)vijk-2]  + 
-Vj+2jkVbi+M  +  8Ut+ljkVb,+1,k  -  SUi-ljkVb^  +  Ui-2 jkVbi.h 


■VkVbj-w 

12A£ 

■Vij+2kVbij+2k  +  Wj+iMij+ik  ~  +  Vv-^vb0.ak  + 

12  At] 


_  •>  ’  'J  ’  - - - - - - - 

“  "  12At? 

-Wijk+ 2vbi]k+2  +  SWjjk+iVbiM  -  8^^-!%,-!  +  Wt]k-2Vb„k- 
~~  12AC 


+ 


( +  v,^j  +  <,^)phit  - 

^ )ijk  +  ^A[-(^.+2jfc  +  Ubi+2jk)wi+2jk  +  8(Ui+ijk  +  Ubt+ljk)wi+ljk 

—8(Ui-ijk  +  Ubi_ljk)roi-ijk  +  {Ui-2jk  +  Ubt_2jk)wi-2jk]  + 

_i— [_(Vij+2fc  +  Vb  )wij+2k  +  8(Vij+ik  +  VbiJ+lk)wij+ik 
12  A  7? 

-8  (Vij-ik  +  Vb^Jwij-n  +  (Vij-2  k  +  Vbijl2k)wij-2k]  + 
12^[-(W,jfc+2  +  Wht]k+2)wijk+2  4-  8(H/,-ji+i  +  Wbijk+1)wijk+ 1 
-8  {Wijk-i  +  Wh^Jwijk- 1  +  {Wijk-2  +  Wbijk_2)wijk- 2]  + 

—  Ui+2jkU>bl+2jk  +  ^Ui+ljkWb^j^  -  8^-ljfc^b,_ijfc  +  Ut-2jk™b^- 

'  12A£ 


+ 


(20) 


0,  (21) 
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- Vij+2kWbiJ±2k  +  %Vij+ikWbii+ik 

12At? 

-WjjwWbijw  +8Wijk+iv>biik+i  -  SW.^-iWb,^!  +  WtJk-iwbijk_2 

‘  “  12AC 


+ 


+ 


<j({'^+’b^+<,^)P>iS‘_(^Al”)ii‘  "  °'  <22) 


Do  the  same  way  as  last  subsection,  we  obtain 


n—  1 


Wiik-4U?ik  +  Viji' 


2A  t 


■  + 


—L-[-Ji+2jk{Ui+2jk  +  Ubt+2jk)Ui+2jk  +  8Ji+ljk{Ui+ljk  +  Ub,+1]k)Ui+l]k 

-L 

—&Ji-ljk{Ui-ljk  +  Ub,_ljk)Ui- life  +  Ji-2jk{Ui-2jk  + 

[-J.j+2fc(Vjj+2fc  +  Vbij+3k)Uij+2k  +  8Jij+lk{Vij+lk  +  Vbij+ik)Uij+lk 

12A?7 

-8J,j_lfc(V1y_u  +  Vb^JUij-lk  +  Jij-2k{Vij-2k  +  Vbti_2k)Ui3-2k\ 

■  ^-[-Jijk+2  {Wtjk+2  +  Wb,]k+2)Uijk+2  +  SJijk+i(Wijk+i  +  Wbljk+1)Uijk+ 1 

—8Jijk-l  {Wijk-l  +  Wbijk_i)Uijk-l  +  Jijk-2{Wijk-2  +  Wbijk_2)Uijk-2] 

-Pi+%sk  +  nri+kik-nPi-h*  +  pi-i» 

+aijk  24  A^ 

1  r  —Ji+2jkUi+2jk  +  16Ji+ljkUi+ljk  —  30 JjjkUjjk  +  16Ji-ljkUi-ljk  ~  Jt-2jk^t-2jk 

12A£2 

-Jjj+2kUjj+2k  ±  Mj+lkUij+lk  -  MJijkUijk  +  16J,-,--ifct/jj-ifc  -  Jrj-TkUjj-M 
~~  —  12At?2 

-Jt7k+2Ullk+2  +  16J,-,-fc+it/fjfc+i  -  30 +  16Jp-fc-iC/,-jfc-i  -  JtJfc-2^^-2i 

ifc' - -  12AC2 

pn-  ps  ,  x  Pj  -  Pb 

f  ScONVijk  (0  +  SviSijk  (0  +  J'..‘Ar+^Ar 

1 

Re 


■Ri{0,iik 

+Pi jk  ‘ 
+7i^‘ 


+  W52^  +  e**53tf*}  =  0  ,  (23) 
3V,'jfc  -  4Vgfc  +  Vjjfc  : 1 


2At 


•  + 


i2Ae 


12At? 


[-Ji+2jk(Ui+2jk  +  Ubi+2jk)Vi+2jk  +  &Ji+ljk(Ui+ljk  +  Ubi+ljk)Vi+1jk 
—8Ji-ljk{Ui-ljk  +  Ub'-ij^Vi-ijk  +  Ji-2jk{Ui-2jk  +  Ub,_2jk)Vi-2jk\ 
[-Jij+2k{Vij+2k  +  Vbtj+2k)Vij+2k  +  8Jij+lk{Vij+lk  +  H,J  +  i  JK'j+U 


-8  Jij-lk{Vij-lk  +  Vbt}_lk)Vl3-lk  +  Jij-2k{Vij-2k  +  Vbi2-2k)Vij- 2k] 

[~Jijk+2(Wijk+2  +  Wb,jk+2)Vijk+2  +  8Jijk+i  {Wijk+i  +  Wbiik+1)V{jk+ i 

-8Jijk-i{Wijk-i  +  WbiiM)Vijk-i  +  Jijk-2(Wijk-2  +  Wbi}k_2)V \jk-2] 


+flijk 


~P.j+|fc  +  27P«j+ifc  ~  27P«i-^  +  P‘7-|fc 
24At? 
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1  -Ji+ijkYi+Tjk  ±  16Ji+iikVi+ijk  -  30 JjjkVijk  ±  lWi-ijkVj-ijk  -  £iz wXtw 

-  Ye{aiik  12Ae2_ 

-Jg+ikVij+Tk  ±  IQJij+ikVij+ik  -  30 JjjkVijk  ±  lQJij-ikVij^ik  -  Jij-2kVij-2k 

12A  If 

—  Jiik+2^ijk+2  +  l6Jijk+lV>jk+l  ~  30J,jfcV|jfc  +  16Jtjfc-lKjfc-l  ~  J x]k-2^'i3k-2  ^ 
+7»'jfc  AC2 


+^l.yfc 


-PF  -  Pw 
A£ 


+  ^20k 


Pj-Pb 

AC 


+  Sccwv, ,*(»/)  +  Svis,Jk  (v) 


-J^{VxhkSlljk+VyijkS2ijk  +  Vzx]kS3,]k}  -  0  ,  (24) 

12Ai 


1 

12  AS 


[-J,+2jfc(t/.+2lfc  +  Ub,+i}kWi+2jk  +  8Ji+ijk{Ui+ljk  +  Ubi+ljkWi+ijk 
-8Ji-ijk{Ui-ijk  +  Ubi_ljk)Wi-ijk  +  Ji-2jk{Ui-2jk  +  Ubi_2jk)Wi-2jk] 


T7T^[-Jii+u{Vij+2k  +  Vbij+2k)Wij+2k  +  8J>j+ik(Vij+ik  +  Vbij+lk)Wij+ik 

l2Arj 

-8Jij-lk{Vij-lk  +  Vbl]_lk)Wij-lk  +  Jij-2k{Vij-2k  +  Vbij_3k)Wij-2k] 

j^[-Jijk+2{Wijk+2  +  Wbtjk+2)Wijk+ 2  +  SJijk+liWiik+i  +  Wbijk+l)Wijk+i 

-SJijk-i{Wijk-i  4-  Wb^^Wijk-i  +  Jijk-2  {W.jk-2  +  Wbijk_2)Wijk-2] 

~Pjjk+  \  +  27PljA.+  i  -  27Pi-k_i  +  Pjjk- 1 
+lijk  24AC 

1  ,  -Ji+2ikWi+2jk  +  lQJi+yjkWi+ijk  -  30 JijkWijk  +  lGJi-ijkWi-ijk  ~  Ji-2jkWi-2jk 
~We^ijk  '  12AS2 

-Jij+2kWjj±2k  +  lQJij+ikWij+1k  -  30 JykWjjk  +  lGJjj-ikWjj-ik  ~  J-.]-2kWtJ-2k 
+ *  ‘  12A772 

-Jijk+2Wijk+2  +  lGJijk+iWjjk+i  ~  30 JijkWijk  ±  lSJijk-iWijk-i  ~  Jijk-2Wijk-2, 

+7  ijk  "  12AC2 

X3~  AS~  +  A2\t  +  Sc0NV,]k^  +  Svis'Jk^ 

+  C,yi}k s2,jk  +c zijkS3ijk)  =  o  ,  (25) 


where 

ScONVijk{4>)  =  J2 -£c{-(Ui+2jk  +  Ubi+2jk)[Ui+2jk{<t>Xijk  -  <j>Xi+2jk) 

JTvi+2jk(4)yl}k  —  fiy.+ijk)  +  wi+2jk{<f>zi]k  ~  0z;+2jifc)] 

+8(Ui+ijk  +  Ubi+1jk)[iii+ijk{<f>xijk  ~  4>xl+ljk ) 

— 8(t/,_ijfc  +  Ubi_l-k)[ui-ijk{4>xijk  ~  01,-iji*) 
~^~‘^i—ljk(^yijk  ~~  ^Vt  —  ijk  )  “I"  Ijki&Zij/c  ^i— 
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+  (U{-2jk  +  U )  [ui-2  jk  —  fixi -2jk) 

+ul-2jA: —  ^Vi-ijk)  Wt-2jk(.<l)Zijk  <Pzt-2jk  )]} 

+  -  1  —  {~(Vij+2fc  +  Vfc0+jk)[Utj+2fe(^iiy*  “  ^xi)+2k) 

12A77 

+t>f j+2fc  {tPtlijk  ~  ^Vij+ik)  wij+2k  [tyztjk  $zi}+?k  )1 

+8(Vij+lJt  +  HtJ+ 1  k  )  [U»  j+ 1 k  (Vxijk  ~  ^Xij+lk) 

~^~Vij+lk{4>Vijk  ~  ^yij+lk)  +  ^ij+lfcC^ZyJr  (t>zi)+ lk)l 

—  8{Vij-lk  +  _  1  fc ) [u«j - 1  * ( fix.jk  ~  ^Xij-ik) 

+vij-\k{4>yXJk  ~  ^Vij-lk)  wij-lk{4>Zijk  ~  ^Ztj-lk  )] 

+  (Vij-2k  +  ) [uij-2k {4>xijk  ~  4>Xij-2k  ) 

+U,-j_2A: _  <f>yn-2k)  +  wv-2k{4>zi]k  ~  <Pz,,_2k)}} 

+  Y^^-{-(W,jfc+2  +  W6ij,+2 ) [Utjfc+2 -  *«iiW)  ■ 

+V*jifc+2(^»i>fc  ~  ^Vijk+2  )  +  wijk+2(4>zijk  ~  tzijk+l)] 
+8(W/,-jfc+i  +  ^6i_,fc+i)[u*jfe+l(^® 

+Vijk+l{<l>yijk  ~  <f>yijk+i)  +  Wijk+lifcijk  ~  faijk+i)] 

-8(W{jk-l  +  Wbijk^i)[uijk-\{4>xi3k  “  4>xijk-i) 

+Vijk-l  (« (t>yijk  ~  <?W- 1 )  +  wijk-l  {4>zijk  “  )] 

+UijJk-2(^Wjfc  “  +  wijk-2{<f>zijk  -  ^o-fc-a)]} 

+  +  +  W^+2>*^i>fc) 

+8t/t’+ljA:('U'6i+ljfc^j:i>?fc  +  Vbi+ijk  ^Vijk  Wbi^ijk<f>zijk) 

—  8£/t_ \jk  iubi_ijk  fixijit  Vbi—ijk  $Vijk  ^bi—ijk^Zijk) 
~\~Ui-2jk{Ubi-2jk(t)Xijk  +  vbi~2jk<f>yijk  "I"  ^bi^jk^Zijk)} 

+  12  A  77  {  ~~  ^+2A:  (^^*J+2fc  ^0*  Vi>ij+2k<l)yijk  Wbij+2k(f)Zijk) 

+8Vij+ik{'U'lij+lk<f>xijk  “I"  Vbij+ik^Vijk  Wbij+ik^Zijk) 
—8Vij-lk{ubij-ik<l)Xijk  +  Vbij-ik^yijk  +  Wbtj^k4>ztjk) 
+Vij-2k{ublJ-2k<t)Xijk  Vbij-ik&Vxjk  ^  Wbij-2k<t>zijk)} 

+  ^2^{  —  Wijk+2{ubijk+2(l>Xijk  "t"  vbijk+2<l)yijk  “t"  Wbijk+Z^zijk) 

+8 Wijk+1  {Ubijk+l  ^Xijk  ^zijk  ) 
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(26) 


SWijk-l{Ubijk-l^ijk  +  Vbijk-\^Vijk  +  Wbijk-l^Zijk) 

+Wijk-2{ubi)k.^xijk  +  vb„k-i4>y,,k  + 

Svisijk (t)  =  ~^^^^~Ui+2jk^k  -  **i+a>fc)  +  Uu'+ljk{4>x'>k  ~ 

+16li,-_l jk(<t>xijk  -  fa.-ljk)  ~  Ui-2]k{<f>xhk  -  <l>Xi-.7ik) 

-Vi+2jk{<f>y,jk  -  <£y;+ 2,fc)  +  16ut+ljfc(0yljfc  -  </>yi+ljk ) 

+  16u,-iji'(<?yl_,k  -  4>y,-llk)  -  u«'-2. jk{4>yijk  ~  Vyi-ijk) 

—  Wi+ljlcifizijk  ~  fizi+ijk)  jki^Zijk  ^Zi+ljk) 

+16Wi-ijk(4>Zijk  ~  fai-lik)  ~  Wi-2jk{4>zi}k  ~  ^*i_a,-fc)] 

+  rJrir~-o[~uij+2k(<t>xijk  ~  <7^+2*)  +  16u,j+ifc(<^xi:,fc  — 
l/A^ 

+  16U{j_lfc(<&Eij*.  —  ^Xij—ik)  ~  ^-tj—2k  {4>x,jk  $Xij—2k) 

-Vij+2k(4>ytjk  -  ^y,'j+2fc )  +  16uij+ifc(<7W  “  ^y.j+ifc) 

+16u,j_ifc(^yt>fc  “  <t>yi}-ik)  *  vij-2k{4>yijk  ~  4>yi,- 2k) 

-Wij+2k(<t>zijk  ~  <f>ziH 2k)  +  16u’>i+i*;(^.^  "  faij+ik) 

+lQWij-ik{<i>zi)k  ~  4>ztJ-lk)  ~  wij-2k{<t>zijk  ~  02i>_2*)] 

+  [— U<i^+2  (^ijk  —  (t>xi]k+2  )  "I”  i^Xijk  <Px%jk+l  ) 

+  16u,jfc-l(<^x,>fc  -  4>Xi,k- 1)  -  Uijk^i^x.jk  -  ^Xijk-l) 

-vijk+2{<j>y,,k  ~  ^ynk+2)  +  16u«'j*+i (^y.jfc  -  ^W+J 
+16Vijk-l{<f>yhk  -  fiyijk- 1)  _  vijk-2{<t>yijk  ~  ijk-i) 

-Wijk+2{<Pzi}k  ~  4>Zijk+  2)  +  16w«jfc+l(&0-*  “  ^.,fc+i) 

+  16iy,jfc-l(^z0jfc  -  <f>Zijk- 1)  ~  Wijk-\{<t>zljk  -  4>zijk- 1  )]}* 

3  Conformal  Mapping 

Since  transition  (especially  in  the  receptivity  regime)  is  very  sensitive  to  any  kind  of  distur¬ 
bance,  physical  or  numerical,  a  highly  accurate  scheme  is  required  for  the  computational  domain 
transformation.  To  assure  the  numerical  scheme  for  the  flow  not  to  lose  accuracy,  we  should  make 
the  accuracy  of  grid  transformation  schemes  at  least  the  same  order  as  that  of  the  flow  scheme 

itself. 

As  the  first  step,  relatively  simple  configurations  are  considered,  and  this  enables  us  to  use  the 
analytical/conformal  mapping  to  obtain  the  required  transformation  coefficients  (Jacobi  metrics). 
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For  an  H-type  grid,  the  Joukowsky  transformation  is  adopted. 

3.1  Conformal  mapping  for  a  quarter  of  cylinder 

For  a  quarter  of  cylinder,  a  three-step  transformation  is  required  to  transfer  the  physical 
domain  to  the  computational  domain  (see  Figure  2): 


1.  (I,??)  =*  (6,t?i): 


t?i  =  V- 


This  will  ensure  the  points  (0,0)  and  (a,  0)  in  the  (£i,f?i)  space  match  the  points  (  a,0) 
and  (0,  a)  in  the  (z,y)  space. 


2.  (Ci,t?i)  =*  (*i»yi): 


Xi 


Vi 


- — — -  —  a, 

aB  +  a-  |fi| 

_ ^CVUnarV} _ 

&CVfnax  "I”  2/1  max  iv  max  m) 


This  includes  stretching  in  both  directions.  The  smaller  the  aB  and  ctc,  the  stronger  the 
stretching  near  the  solid  wall  and  the  leading  edge,  respectively.  The  associated  inverse 
transformation  is 


TJmaxVl  {PC  d"  yi  max  ) 

VlmaA^C  +  J/l) 


3.  (x?y)*  Joukowsky  transformation 


x  a _ 

Xl  ”  2  +  x2  +  y2 

y  a2 

yx  ~  2  ^  x2  +  y2 


The  associated  derivatives  in  this  step  are 

1  a2  y2  —  x 2 
Xl*  =  2  +  T(:r2  +  y2)2’ 

a2xy 

Xl»  =  _(*2  +  y2)2’ 
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(28) 


=  g-4a£^g  +  ^I[glyy((7B  +  a;i  +  o)  -2»iy]. 


•  Xi  <  —a: 


c 

sx 


£y 

?xx 

£yy 


( OB  +  ^1  +  a) 


Q_A^Bi^B_±a) 

(ctb  ~  (X1  +  a))2 
oao_b{qb_  +  a) 

(<7b  ~  (^1  +  a))2 

+  °L{  {aB  _  (I,  +  „)]  +  to?.}, 

{ob-  (*i  +  «)r 

WBWB  +  *} Ux  [oB  _  (*,  +  „)]  +  2xl  }. 

{°B  -  (*i  +  a)) 


\2XlV’ 


(29) 


normal  direction: 


rirnaxjvc  ±  yirnax)  ?cy  lx 
Vx  ~  yimax  (^C+yi)2’ 

_  Vmaxi^C  d"  yimax ) 

Vy  "  yimax  (^c  +  yi)2’ 

v  =  ?£Vmax(°C_  +  Kwly^  (gp  +  ;/lmo  J  -  2y\x (eC  +  Vl? , 

y  i  mfli 

nw  =  ac7?mal(ac  +  yi~— yiy>c  +  yimax)  -  2y?v(frc  +  yi)3- 

2/lmax 

For  the  geometry  other  than  a  cylinder,  more  transformation  steps  are  required. 
3.2  Conformal  mapping  for  the  leading-edge  flat-plate  juncture 
The  transformation  process  can  be  described  as  follows  (see  Figure  3): 


(30) 


i.  (£,  y)  =$”  (£o>yo): 


£o  =  —  t 

°A 

■no  =  y  ■ 


2.  (£o,  Vo)  ^  (£i>yi)- 


6 

m 


Qb£o  _ 

OB  +  a  ~  ICol 

_ Q£V  ImaxVO _ 

OCVmax  +  V  1  max  (y  max  yo) 
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V 

% 

m 

y  i 

y 

Tim  n  r 

Vmax 

v/l  max 

Umax 

Vlmax 

a 

ii 

r 

3  _ 

- -  7 

^0 

— 

a  -cl  ( 

)  -a  U 

—  «:ro  u 

Figure  3.  Transformation  process  for  an  elliptic  leading-edge. 


3-  (£i,t?i)  (*i,yi):  Joukowsky  transformation 


6 

Vi 


^l  +  yj 


n 

a * 

*5  +  y? 


). 

)• 


4.  (xi,t/i,zi)  =J>  (z,y,z): 


x  = 

y  = 


CTDXl[l  + 


fc2  -  1  1  . 

(fc  +  l)2  x2  +  y2 
fc2  -  1  1  , 

(k  +  l)2  x\  +  yf 


where  axes  of  the  ellipse  are  (ro,  &ro),  with 

ro{k  +  1) 

a  =  — - - 

2&D 


3.3  Conformal  mapping  for  Joukowsky  airfoil 


The  transformation  process  for  Joukowsky  airfoil  can  be  described  as  follows: 

1.  (€,v)  =*  (£oi  Vo): 
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2.  (Co,  Vo)  =k  (6>  7?i): 


6  = 


°b£, o 


-  a, 


GB+a~  ICo| 

_  _ gCgimox^O _ 

CFcVmax  +  VlmaxiVmax  ~  Vo) 


3-  (Cii7?!)  =>■  (xi,yi):  Joukowsky  transformation 


Cl  =  T-(l  + 


2'"  *i  +  tfi 

j/i  a2 

711  ~  2  1  x?  +  y\ 


), 


)• 


4.  (xi,t/i)  =>■  (X2,  2/2) : 


kfa 

2:2  =  Xl"F+T’ 

2/2  =  2/i- 


5.  (x2,  j/2)  -*  (z,y): 

,  a2  1  x 

1  =  °EXl{1  +  W+WW+y!  ’ 

,  °2  1  x 

y  =  <TE!/2(1-  (fc'+l)=xl  +  y|)- 

Here,  <  1  is  used  to  change  the  ratio  of  half-thickness/chord-length  of  the  airfoil.  For  example, 
with  k'  =  0.0116,  the  ratio=l/16.  Also,  we  use  oe  to  normalize  the  airfoil  and  make  the  half¬ 
thickness  equal  to  1  (see  Figure  4). 

The  chain-rule  is  applied  to  obtain  the  required  metrics.  For  example, 

ix  =  Cn^i*  +Ct/i2/ix  +Czi2ix> 

Cy  =  Cxi^ly  +  CyiJ/ly  +  Czi^ly) 

Cz  =  Cxi^i*  +  Cyi  2/iz  +Czi  Zizi 
T)x  =  »?*1*lx  +  »?yi  yi«  +7?*i2lz) 

JJj,  =  7?XlZls,  +  7?yiyiv  +7?z12lS),  . 

^  =  Vx^iz  +  TlyiVlz+Vz^u, 

Cx  =  Cxi  XU+  Cyi  yi*  +  Czj  Z\x ) 

Cy  =  Cxi*ly  “k  Cyi I/ly  ~k  Czi^ly) 

Cz  =  CxiXl2  4"  Cyi  2/lx  "k  Czi  ^12) 
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Figure  4.  Transformation  process  for  a  Joukowsky  airfoil. 

with 


f*l 

= 

ZtlZlXl  +  +  CClClXi? 

€yi 

= 

C^iCiyi  “^CCiCiyn 

Cm 

= 

"1“  Cm7?!*!  “^CCiCl^ii 

Vxi 

= 

^iCixi  +  VmVixi  +  %Cixn 

Vyi 

=: 

^lClyi  “t“  VmVlyi  ^CiClyi > 

*1*1 

= 

+  ^m7?^  +  %iCl*n 

Cx  1 

= 

C^i^ixi  “t“  Cm7?!*!  “1”  CciClxn 

Cvi 

= 

Ce^iyi  +Cm7?iyi  +  CciCiyn 

Czi 

= 

CciCui  +Cm7?i2i  +CciCui* 

The  second-order  derivatives  can  be  obtained  through: 

£xx  =  £x£^x  +  £xt)Vx  ~h  £*<Cxt 

£yy  —  +  ZyrtVy  +  ^vcCvi 

fizz  —  4  CzTJ^Z  4  Cz(Czi 

=  Vx^x  4"  VxrtVx  4  VxC.Qxi 
Tjyy  -  Vy^y  +  VyriVy  +  VyiCyi 
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Vzz  =  Vziiz  +  VzriVz  +  VztCz, 
Cxx  —  Cx{£t  +  CxvVx  +  CxcCx, 
Cyy  =  Cyt£y  +  CynVy  +  CycCjn 
C zz  —  C  ziiz  +  (zyVz  +  GcCz- 


4  Numerical  Process 

In  the  computational  (£,  t?,  C)  space,  the  grids  are  uniform.  Suppose  U,  V,  W  are  defined  in 
terms  of  a  staggered  grid  in  the  computational  space  (see  Figure  5).  Here,  the  values  of  u,  v,  w,  P 
are  associated  with  its  cell  centers,  U  with  centers  of  the  cell  surfaces  parallel  to  the  (t?,C)  P^ne, 
V  with  centers  of  the  cell  surfaces  parallel  to  the  (£,C)  plane,  and  W  with  centers  of  the  cell 

surfaces  parallel  to  the  (£,77)  plane. 


Figure  5.  Staggered  grid  structure  in  the  computational  (£,??,  C)  space. 

The  approaches  used  in  this  work  include 
•  Base  flow 

1.  Hybrid  scheme  in  space  (Patankar  &  Spalding,  1972) 

2.  Line-distributive  relaxation  (Liu  &;  Liu,  1993) 
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3.  Multigrid  FAS  (Brandt,  1984) 


•  Perturbation  flow 

1.  Second-order  backward  Euler  in  time 

2.  Fourth-order  central  difference  in  space  (Liu  k  Liu,  1993) 

3.  Line-distributive  relaxation  (Liu  k  Liu,  1993) 

4.  Buffered-outflow  boundary  condition  (Liu  k  Liu,  1993) 

5.  Inter-grid  dissipation  (Liu  k  Liu,  1995) 

A  basic  state  flow  field  is  first  solved.  Since  the  H-type  grid  is  used,  we  do  not  need  large 
number  of  grid  points  in  the  wall-normal  direction  (C-type  grid  usually  needs  more  grid  points 
because  the  wall-normal  direction  is  almost  parallel  to  the  main  flow  direction  in  some  region) . 
The  perturbation  system  is  then  solved. 

5  Computational  Results 

Based  on  the  above  formulation,  a  number  of  computer  codes  for  grid  generation,  base  flow, 
and  perturbation  flow  have  been  developed. 

5.1  Code  validation 

5.1.1  Linear  stage 

The  validation  of  the  code  is  first  examined  on  the  linear  instability  stage  for  a  smooth  flat 
plate. 

The  inflow  perturbation  for  the  2D  case  is 

u(0,y,  z,t)  =  Real{e2<f</>2^(!/)e  ,W2dt},  (31) 

for  the  3D  case  is 

„(<>,!,, z,t)  =  <32) 

and  for  the  3D  nonlinear  breakdown  case  is 

u(0,  y,  z,  1)  =  Real{<M4‘>  (y)^'  +  ^4’ (33) 
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Here,  e2d  and  f3d  are  the  amplitude  of  2D  and  3D  disturbances.  The  grid  transformation  function 
in  the  wall-normal  direction  is  given  by 

2,(7?)  = - -JtZZEl - (34) 

T)max&  ”1”  Umaxylmax  Vj 

Here,  we  choose  a  =  3.5,  and  ymax  =  Vmax  =  50  (based  on  8  ). 

Both  2D  (0  =  0)  and  3D  (0  ^  0)  cases  are  studied.  Typical  3D  results  are  shown  in  Figures  6 
y  ^  17Q  x  34  x  18  grid  is  employed.  We  divide  the  whole  computational  domain  into  a  physical 
domain  and  a  buffer  domain,  the  latter  is  one  T-S  wavelength  long.  Let  Re *  =  688,  u2d  =  u3d  = 
u  =  0.075 {F  =  uj/Rer  x  106  =  109),  and  0  =  0.2.  The  comparison  of  the  amplitude  distributions 
in  both  the  wall-normal  direction  (eigenfunctions)  and  streamwise  direction  (amplifications)  shows 
a  good  agreement  with  the  LST  results. 


CASE 

Fundamental  (Figure  7  (a)) 

Re m 

1221.77 

grids 

290  x  34  x  30 

length 

16  \ts  physical  domain,  2  Ajs  buffer 

height  (based  on  6*) 

50 

time-step/T-S  period 

220 

inflow  disturbance 

e2d  =  0.0056,  e3d  =  0.0004 

other  parameters 

0  =  0.2451 
aR  =  0.2451 
aR2d  =  0-2690 
u2d  =  <+>3d  -  0.0928 
cr  =  4.05 

CASE 

Subharmonic  (Figure  7  (b)) 

Re * 

1221.77 

grids 

338  x  34  x  26 

length 

19  A ts  physical  domain,  2  A xs  buffer 

height  (based  on  5*) 

50 

time-step/T-S  period 

220 

inflow  disturbance 

e2d  =  0.003,  e3d  =  0.000125 

other  parameters 

0  =  0.178 

Oiji  =  0.12965 
aR,d  =  0.2690 

U2d  =  0.0928,  u3d  =  0.0464 
a  =  4.15 

Table  1.  Parameters  for  transition  over  flat  plate. 


5.1,2  Whole  process  of  transition  over  a  smooth  flat  plate 

To  verify  the  ability  of  this  code,  further  tests  on  a  smooth  flat  plate  were  performed.  Using 
the  experimental  input  data  of  Saric  et  al.(1984),  we  did  simulations  for  both  the  fundamental 
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breakdown  and  subharmonic  breakdown.  The  inflow  boundary  condition  is  assumed  to  have  the 
form  of  formula  (33).  The  computational  parameters  are  summarized  in  Table  1,  and  the  results 
are  shown  in  Figure  8.  The  numerical  results  match  the  experimental  results  very  well. 

5.2  Receptivity  of  flat-plate  with  an  elliptic  leading-edge 

In  this  case,  we  set  the  Reynolds  number  based  on  the  small  axis  of  the  ellipse  Ren  =  2400. 
Figure  9  depicts  the  configuration  considered  here. 

The  base  flow  is  characterized  by  large  streamwise  and  normalwise  velocity  gradients  around 
the  leading-edge,  and  appearance  of  a  maximum  adverse  pressure  gradient  exactly  at  the  leading- 
edge-flat  plate  juncture.  The  boundary  layer  adjusts  to  recover  to  freestream  pressure  after  an 
acceleration  over  much  of  the  nose  surface  (mainly  due  to  a  strong  favorable  pressure  gradient 
there).  Figure  10  depicts  results  of  the  current  base  flow  pressure  gradient  along  the  stagnation 
line  and  the  wall  surface,  and  that  of  Buter  and  Reed  (1994)  for  comparison.  They  are  in  good 
agreement.  The  adverse  peak  value  of  the  surface  pressure  gradient  caused  by  body  curvature 
changes  at  the  juncture  is  considered  to  be  one  of  the  receptivity  sources  other  than  leading-edge 
(Goldstein  k  Hultgren,1989).  It  should  be  noted  that  a  very  small  but  almost  constant  adverse 
pressure  gradient  will  amplify  the  instable  mode  once  generated  and  after  reaching  the  first  branch 
of  its  neutral  curve.  Normal  components  of  the  base  flow  at  the  juncture  and  exactly  the  same 
downstream  station  shown  in  Figure  11  are  found  to  agree  well  with  that  of  Buter  k  Reed  (1994). 
The  boundary  layer  exhibits  clearly  a  non-parallel  characteristic  which  is  not  accounted  for  in  the 
classic  linear  stability  analysis. 

The  unsteady  flow  resembles  a  temporal  and  spatial  evolution  after  an  external  harmonic 
vortical  disturbance  is  imposed  at  the  upstream  boundary, 

u{x0,y,t)  =  Real{een/orceA(|y|  -  yc)e~ (MVc)  e~lujt },  (35) 

where  A  is  a  parameter  for  adjusting  the  amplitude  of  enforced  disturbance,  u  is  the  frequency  of 
the  imposed  disturbance,  k  is  used  to  change  the  shape  of  u’s  profile,  and  yc  is  the  central  point 
of  the  disturbance.  Here,  we  use 

u  =  1.0,  k  =  0.0075,  t/c  =  0.15. 

Note  the  disturbance  is  bounded  very  closely  to  the  stagnation  line  of  the  base  flow,  as  T-S 
waves  are  expected  to  be  generated  by  a  source  term  that  is  roughly  the  product  of  a  Stokes 


23 


shear  layer  adjacent  to  the  wall  and  the  streamwise  gradient  of  the  mean-flow  velocity  (Goldstein 
&  Hultgren,  1989).  The  boundary  layer  does  involve  a  triple-deck  structure  after  passing  the 
leading-edge-flat  plate  juncture,  and  this  can  be  seen  from  the  instantaneous  streamwise  velocity 
component  distributions  at  different  distances  from  the  wall  after  the  10th  cycle  of  forcing  (Figure 
12).  The  response  is  evidently  a  complex  combination  of  signals  of  varying  wavelengths.  The 
nature  of  the  flow  near  the  wall  is  influenced  almost  exclusively  by  the  presence  of  T-S  waves, 
while  in  the  main  inviscid  region  of  the  boundary  layer  convected  disturbances  dominate.  Figure 
13  shows  the  spatial  evolution  of  the  generated  T-S  wave  (u  component)  aft  of  the  pressure 
gradient  maximum(at  the  juncture).  The  computed  location  of  Branch  I  is  about  x  =  11.0.  A 
comparison  of  T-S  amplification  factor, 

N  =  log(A/Ao) 

along  the  flat  plate  with  that  obtained  by  Buter  k  Reed  (1994)  is  depicted  in  Figure  14.  The 
current  result  is  not  as  smooth  as  Buter  k  Reed,  and  one  possibility  is  our  input  disturbance  is 
not  exactly  the  same  as  theirs.  Also,  from  our  simulation,  we  find  the  original  convection  wave 
(long  wavelength  wave)  damps  not  as  fast  as  that  of  Buter  k  Reed.  Since  this  wave  has  the  same 
frequency  as  the  excited  T-S  wave,  the  amplification  factor  has  some  oscillations  which  involves 

the  effects  of  the  damping  long  wave. 

The  ratio  of  generated  T-S  wavelength  to  that  of  freestream  disturbance  wavelength  is  found 
to  be  around  0.40,  and  quite  agreeable  to  the  estimation  by  Goldstein (1983),  which  is  0(e)  where 
e  =  0.39  if  calibrated  by  Re~hl/&.  The  receptivity  coefficient  is  defined  here  as  a  convention  by 
adopting  the  maximum  rms-value  of  the  disturbance  at  a  location  far  upstream  of  the  boundary 
layer.  The  level  of  receptivity  in  the  current  study  at  Branch  I  is  2.05%,  while  the  same  parameter 
in  Buter  k  Reed  (1994)  is  reported  to  be  around  2.86-3.30%.  This  discrepancy  may  be  due 
to  the  quantification  of  receptivity  levels  in  the  manner  adopted,  as  the  level  of  receptivity  is 
more  inherently  related  to  some  freestream  disturbance  which  is  varying  dramatically  along  the 
boundary  layer  (Buter  k  Reed, 1994).  Note  that  the  receptivity  coefficient  of  the  same  geometry 
to  freestream  sound  in  Lin  et  al.{ 1992)  is  of  the  order  10-1,  much  higher  than  both  the  above  two 
computations. 

5.3  Receptivity  of  2D  and  3D  Joukowsky  airfoils 

The  upstream  boundary  condition  of  the  base  flow  at  infinity  is  considered  to  be  with  zero 
attack  angle,  and  the  external  freestream  disturbance  is  introduced  symmetrically  to  the  axis  of 
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the  Joukowsky  airfoil.  This  enables  us  to  carry  out  simulations  within  the  upper  half  of  the  flow 
field.  Actually,  all  2D  and  3D  results  are  obtained  by  the  same  code  by  assuming  periodic  spanwise 
boundary  condition.  The  ratio  of  the  half-thickness(fi)  of  the  airfoil  to  its  chord-length(I)  is  1/40, 
with  the  crest  location  about  1/4  of  the  chord-length.  Other  parameters  for  2D  and  3D  simulations 

are  listed  in  Table  2. 

Similar  to  that  of  the  flat  plate  with  an  elliptic  leading-edge,  the  base  flow  of  the  present  case 
has  a  strong  favorable  pressure  gradient  near  the  stagnation  point;  the  boundary  layer  remains 
very  thin  until  reaching  the  adverse  pressure  gradient  regime.  Except  for  a  striking  distinction 
near  the  trailing-edge  of  the  airfoil,  of  which  the  influence  we  shall  discuss  later  on,  the  present 
base  flow  differs  from  the  former  mainly  in  the  surface  pressure  gradient  distribution  along  the 
wall.  Shown  in  Figure  15,  the  surface  pressure  gradient  is  positive  everywhere,  with  its  maximum 
value  rather  smaller  than  that  in  Figure  10.  This  is  because  the  leading-edge  aspect  ratio  is 
smaller  (recall  a  comparison  of  1/10  with  1/6),  and  the  surface  curvature  discontinuity  has  been 
removed  by  smoothness  of  the  Joukowsky  airfoil. 


2D  case 

3D  case 

Rei 

1.6  x  10s 

1.8  x  10b 

L/h 

40 

36 

grids 

482  x  50  x  4 

402  x  32  x  26 

domain  height 

25  h 

25  h 

time-step/cycle  of  force 

200 

200 

disturbance  amplitude 

^enforce  =  0.01 

C2d  —  0.02, 

63d  =  0.02 

frequency 

u=  1.0 

u  —  1.2 

(F  =  250) 

(F  =  240) 

Table  2.  Parameters  for  receptivity  of  2D  and  3D  Joukowsky  airfoils. 

It  is  reasonable  to  guess,  in  this  case,  T-S  waves  will  be  present  later  along  the  wall  and 
undergo  a  longer  period  of  time  of  attenuation  before  enlarged.  Figure  16(a)  demonstrates  the 
process  of  this  development  along  a  fixed  grid  line,  j  =  5,  roughly  equivalent  to  a  height  of  0.08 h 
above  the  surface.  T-S  waves  become  observable  only  in  the  rear  half  of  the  airfoil,  in  contrast  to 
an  immediate  appearance  aft  of  the  leading-edge  juncture  in  the  flat  plate  case.  From  Figure  16(b) 
we  find  the  same  complex  combination  of  varying  wavelength  through  the  airfoil  boundary  layer, 
and  the  wavelength  conversion  as  well.  Again,  the  ratio  of  T-S  wavelength  to  that  of  freestream 
disturbance  is  found  to  be  about  0.40-not  claimed  here  accurately  because  of  variations  of  T-S 
wavelength  itself. 
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With  deferred  onset  of  amplification  of  the  most  instable  wave,  the  perturbation  flow  has  a 
N  factor  (Figure  17)  relatively  small  along  the  body  surface.  It  grows  much  fast  thereafter,  and 
keeps  increasing  throughout  the  wake.  This  coincides  with  the  above  analysis  of  the  disturbance 

velocity. 

The  3D  simulation  is  complemented  by  introduction  of  an  external  spanwise  spatial  harmonic 
disturbance(see  Table  2)  at  the  upstream  boundary.  From  Figures  18  and  19  we  can  find  that  not 
only  the  2D  waves  can  be  recepted  after  the  disturbances  pass  the  leading-edge  region,  but  also 
the  3D  oblique  modes.  The  difference  is  that  there  is  no  wavelength  conversion  in  the  spanwise 

direction. 

5.4  The  whole  process  of  transition  around  3D  Joukowsky  airfoils 

In  this  case,  a  Joukowsky  airfoil  with  infinite  width  and  zero  attack  angle  is  considered.  The 
spanwise  boundary  condition  is  thus  assumed  to  be  periodic. 

5.4.1  Enforced  transition 

We  first  investigate  the  process  of  enforced  transition  in  a  2D  Joukowsky  airfoil  boundary 
layer.  In  this  case,  the  disturbance  is  imposed  into  the  boundary  layer  directly.  The  associated 
parameters  for  numerical  simulation  are  shown  in  Table  3,  and  the  grids  in  (x,  y)-plane  are  shown 

in  Figure  20. 

The  simulation  was  carried  out  on  an  IBM  RS  6000/590  workstation,  and  required  approxi¬ 
mately  100  hours  of  CPU  time  and  65  MB  of  memory. 


ReL 

126,000(^  =  3,500) 

L/h 

36 

grids 

298  x  34  x  26 

domain  height 

25  h 

time-step/forcing  cycle 

250 

disturbance  amplitude 

e7d  =  0.01,  e3d±  =  0.005 

other  parameters 

u  =  0.7(F  =  228.6  x  10~6) 

(3=  1.0 

location  of  forcing 

x / h  —  8.19  (from  the  leading-edge) 

points  on  the  airfoil 

214  (stream wise) 

Table  3.  Parameters  for  simulation  of  the  whole  process  of 
transition  over  a  Joukowsky  airfoil  (enforced  transition). 


The  periodic  disturbance  is  imposed  at  xjh  =  8.19  (measured  from  the  leading-edge).  Though 
we  still  use  formula  (33),  the  distribution  of  the  perturbation  in  the  wall-normal  direction  can  no 
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longer  be  called  “eigenfunction” .  Figure  21  depicts  the  instantaneous  contour  plots  of  perturbation 
ux,  uz,  and  |u|  on  the  j  =  6  (»  0.07A  from  the  wing  surface)  grid  surface  after  the  10th  cycle  of 

forcing. 

Figure  22  gives  a  response  to  the  flow  field  which  is  subject  to  the  periodic  disturbance  input 
at  the  imposed  location.  It  is  quite  clear  that  as  the  perturbation  waves  travel  downstream,  the 
nonlinear  affection  increases.  At  some  downstream  location,  the  sign  of  randomness  appears,  and 
the  so-called  spikes  can  be  observed.  The  flow  field  is  no  longer  periodic. 

5.4.2  Natural  transition 

In  this  case,  a  similar  configuration  as  in  the  above  subsection  is  investigated.  Now,  the 
disturbance  is  not  imposed  directly  into  the  boundary  layer,  but  from  free-stream.  The  pattern 
of  disturbance  is  still  the  same:  one  2D  wave  and  two  oblique  waves.  The  parameters  used  in  this 
simulation  are  listed  in  Table  4. 


ReL 

200,  000(Re^  =  5000) 

L/h 

40 

grids 

402  x  32  x  26 

domain  height  (based  on  h) 

25 

time-step/cycle  of  forcing 

200 

disturbance  amplitude 

e2d  =  0.02,  e3d  =  0.02 

other  parameters 

u  =  1.2(F  =  240  x  10"6) 

/ 3=  1.0 

points  on  the  airfoil 

322  (streamwise) 

Table  4.  Parameters  for  simulation  of  the  whole  process  of 
transition  over  a  Joukowsky  airfoil  (natural  transition). 


The  simulation  was  carried  out  on  an  IBM  RS  6000/590  workstation,  and  required  approxi 
mately  150  hours  of  CPU  time  and  124  MB  of  memory  to  run  18  forcing  periods. 

The  periodic  disturbance  is  imposed  at  xjh  ss  — 15  from  the  leading-edge.  Figure  23  depicts- 
the  instantaneous  contour  plots  of  perturbation  ux,  uz,  and  |w|  on  the  j  =  5  («  0.07/r  from  the 
wing  surface)  grid  surface  after  the  17th  cycle  of  forcing.  Figure  24  gives  a  response  to  the  flow 
field  which  is  subject  to  the  periodic  disturbance  input  at  inflow.  It  is  quite  clear  that  different 
from  the  previous  case,  this  case  includes  a  wavelength  conversion  process.  According  to  Buter 
and  Reed  (1994),  the  receptivity  coefficient  for  the  freestream  vorticity  receptivity  is  affected  by 
the  aspect  ratio  of  the  leading-edge  and  continuity  of  surface  curvature.  In  this  study,  there  is  no 
curvature  discontinuity  on  the  wall  surface,  and  the  aspect  ratio  of  the  leading-edge  is  quite  large 
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(about  1:10),  thus  the  receptivity  coefficient  is  very  small.  Nevertheless,  we  still  observed  the 
nonlinear  breakdown  at  some  downstream  location:  the  sign  of  randomness  appears,  the  so-called 
spikes  can  then  be  observed,  the  flow  field  is  no  longer  periodic  in  time,  and  the  vortex  breakdown 

appears. 

Figure  25  depicts  the  contour  plots  of  total  u)z  at  both  peak  and  valley  planes  during  the  17th 
forcing  period.  Though  the  process  is  quite  similar  to  that  of  the  previous  case,  the  breakdown 
location  was  delayed  because  the  amplitude  of  effective  T-S  wave  induced  from  receptivity  is  very 

weak. 

There  is  no  doubt  the  current  numerical  simulation  has  already  provided  a  clear  picture  of  the 
transition  process  with  relatively  coarse  grid  and  medium  Reynolds  number. 
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Figure  7.  (a)  Comparison  of  the  maximum  u,  u,  and  w  amplitudes  with  LST.  (b)  Comparison  of 
the  numerical  and  LST  profiles  of  velocity  at  x/6m  =  114.91  (3D  case). 


2 

Figure  8.  (a)  Comparison  of  normalized  perturbation  velocity  amplitude  against  Res{=  Rex  ) 
with  experimental  results  of  Saric  et  al.  (1984)  for  the  fundamental  breakdown  case,  (a)  Com¬ 
parison  of  normalized  perturbation  velocity  amplitude  against  Res  at  *  =  0.4  with  experimental 
results  of  Saric  et  al.  (1984)  for  the  subharmonic  breakdown  case. 
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Figure  9.  Sketch  of  a  flat  plate  with  an  elliptic  leading-edge. 


Figure  10.  Base  flow  surface  pressure  gradient  of  the  flat-plat  with  an  elliptic  leading-edg< 


Figure  11.  Base  flow  normal  velocity  profile  at  x=0,  and  x=19.2. 
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12.  Instantaneous  perturbation  u  distribution  after  10th  cycle  of  forcing 
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Figure  13.  Spatial  evolution  of  generated  T  -  S  waves. 


current  results 


*  ♦  LBT  results  (from  Buter  k  Reed,  1994) 

*  *  DN0  results  (Buter  k  Reed,  1994) 
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Figure  14.  Computed  perturbation  amplification  factor. 
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Figure  15.  Base  flow  surface  pressure  gradient  for  the  2D  Joukowsky  airfoil. 


Figure  16.  (a)  Spatial  evolution  of  generated  unsteady  waves  along  fixed  grid  line  after  10th 
cycle  of  forcing,  (b)  Instantaneous  perturbation  u  distribution  along  two  grid  lines  after  14th 

cycle  of  forcing. 


Figure  17.  Computed  perturbation  amplification  factor  for  the  2D  Joukowsky  airfoil. 
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Figure  20.  Grids  for  the  Joukowsky  airfoil  in  the  (x,  y)-plane. 
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Figure  21.  Instantaneous  contour  plots  of  perturbation  ux,  uz,  and  |u;|  on  the  j  —  6  (~  0.06 h 
from  wall  surface)  grid  surface  at  t  =  10T  for  the  enforced  transition  case. 
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Figure  24.  Instantaneous  disturbance  velocity  components  at  t  =  17T  along  a  grid  line  (~  0.7 h 
from  the  wing  surface)  for  the  natural  transition  case. 


Total  uz  on.  the  z  =  — 1.701  /  plane  at  t  —  16T,  16.5T,  17T. 


Figure  25.  Instantaneous  contour  plots  of  total  uz  on  the  z  —  0.1309  and  z  —  1.701  f 

planes  during  the  17th  forcing  period  for  the  natural  transition  case. 


6  Concluding  Remarks 

From  the  above  discussion,  we  can  conclude 

•  The  current  numerical  approach  can  handle  fairly  complicated  geometries  by  using  the  high- 
order  grid  transformation. 

•  The  wavelength  conversion  has  been  numerically  simulated,  showing  out  capability  to  sim¬ 
ulate  more  practical  cases  of  natural  and  enforced  transition. 

•  3-D  simulation  for  receptivity  and  transition  basically  works.  We  can  use  a  rather  coarse 
grid  to  study  the  behavior  of  most  energetic  eddies,  and  describe  some  physics  of  transition 
which  may  help  the  experimental  and  theoretical  works. 

•  For  the  realistic  problem,  the  memory  size  and  CPU  cost  will  increase  further,  so  a  code 
under  parallel  architecture  needs  to  be  developed  soon. 
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